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Abstract 

Gravitational perturbations which are present in any realistic stellar collapse to a black hole, die off in the exterior of 
the hole, but experience an infinite blueshift in the interior. This is believed to lead to a slowly contracting lightlike 
scalar curvature singularity, characterized by a divergence of the hole's (quasi-local) mass function along the inner 
horizon. 

The region near the inner horizon is described to great accuracy by a plane wave spacetime. While Einstein's 
equations for this metric are still too complicated to be solved in closed form it is relatively simple to integrate them 
numerically. 

We find for generic regular initial data the predicted mass inflation type null singularity, rather than a spacelikc 
singularity. It thus seems that mass inflation indeed represents a generic self-consistent picture of the black hole 
interior. 



1 Introduction 

The strong cosmic censorship hypothesis suggests the ex- 
istence of a central spacelike curvature singularity at the 
core of the black hole. The Kerr solution which describes 
the exterior of the hole very accurately however possesses 
a timelike singularity in its interior fjj. This singularity 
lies behind an interior Cauchy horizon, which in the Kerr 
case coincides with the inner horizon of the black hole. 
Penrose argued on rather general grounds that this inner 
horizon should be unstable Perturbations, present 
in any realistic stellar collapse S, falling into the black 
hole experience an infinite blueshift as they near the inner 
horizon. At the same time scattered radiation crossing 
the Cauchy horizon cause its contraction. These two ef- 
fects lead to the formation of a lightlike scalar curvature 
singularity along the Cauchy horizon ^| which slowly 
tapers on to the final spacelike central singularity. While 
the latter is believed to be of the oscillatory BKL type || 
the former is rather weak and smooth in the sense that is 
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fully characterized by the divergence of a single function 
along the Cauchy horizon, which can be identified with 
the (quasi-local) mass function of the black hole. Hence 
one speaks also of a mass inflation singularity. 

Mass inflation is well established in the spherical case, 
both analytically ^, || and numerically @ . The analysis 
of the generic, non-spherical case on the other hand is 
much more difficult. Earlier attempts were either some- 
what inconclusive Q or assumed additional properties 
like slow rotation || . Nevertheless these analyses all seem 
to agree with the spherical mass inflation picture. 

Only recently has the full non-spherical problem been 
tackled [[To] . In this latest approach the full set of Ein- 
stein's equations are expanded in an inverse curvature 
series, utilizing a recently developed double null formal- 
ism 0. 

However the equations governing the shear (which en- 
codes the gravitational perturbations leading to mass in- 
flation) are, not surprisingly, a pair of nonlinear coupled 
partial differential equations, even to lowest order in the 
above mentioned expansion. This makes the analysis of 
the structure of the Cauchy horizon singularity exceed- 
ingly difficult. One way to handle this problem is a nu- 



merical analysis of the equations. 

The purpose of this paper is to report the results 
of such a numerical experiment. The numerical solu- 
tion of the Einstein's equations, supplied with the correct 
boundary conditions indeed confirms the analytic predic- 
tion and thus places the mass inflation scenario on firm 
ground. 

2 The plane wave model 

The investigation of the Cauchy horizon is simplified by 
two facts. First, descent into a black hole is necessarily 
progress in time. Remember that the radial coordinate r 
in the Kerr case becomes timelike in the interior. Thus 
we don't need to know the (quantum) physics of central 
core to analyze the outer layers of the black hole interior. 




Figure 1. The Penrose diagram of a stellar col- 
lapse leading to a black hole. The region above 
the potential barrier (gray) is modeled to a high 
accuracy by a plane wave space time. The dashed 
lines represent the characteristic initial surfaces 
used to numerically solve Einstein's equations. 
Note that the "point" H is not part of the space- 
time, but is an artifact of the mapping. 

Second the high blueshift of the infalling radiation 
only occurs behind the interior gravitational potential 
which scatters this radiation once again §]. Thus in- 
falling radiation will propagate on a nearly stationary 
background. Only after the scattering occurs will the 
blueshift set in at which point the back reaction onto 
the spacetime can no longer be neglected. This leaves 



us effectively with a cross flow situation near the Cauchy 
horizon. Figure 1 summarizes the situation. 

Let us now focus our attention to this inner region 
near the Cauchy horizon. The most general metric we 
can write down is of the form [fill 

ds 2 = ~2e x dudv + g ab (d6 a + s a du) (d6 b + s b du) (1) 

where the lower case indices run from 2 to 3 and 9 a denote 
spacelike coordinates. A more detailed discussion of a 
metric of the form ([!]) appears in |ll[] . 

To get some idea of the expected behaviour of the 
metric coefficients in (Q) let us draw the comparison to 
the spherical case described in more detail in ||. For a 
spherical line element we have s a = and g a bd0 a d9 b = 
r 2 d£l 2 in the usual polar coordinates. We therefore would 
expect g a b to stay finite as the Cauchy horizon is ap- 
proached. The important point to note is that r stays 
finite at the Cauchy horizon. Thus we do not deal with 
a r = type singularity of the type expected deeper in 
the hole. 

The gravitational perturbations of the Reissner-Nord- 
strom spacetime which, because of sphericity, must be 
modeled by a scalar field, decay as an inverse power law 
in the advanced and retarded coordinates u and v [^|. 
In these coordinates we have exp(— A) ~ expKo(w + v )- 
Raychauduri's equation then dictates 

r. v ~ v~ q and r „ ~ v~ p , (2) 

so that 

27* iV 7\ u e — g^ 1 t^t v 

= 1 - 2m{ ^ v) + J ~ ft- 

as v — > oo. Here kq is the surface gravity of the unper- 
turbed inner horizon and p = q — 2 is related to the multi- 
pole moment I of the infalling perturbation by p = Al + 4. 

This divergence signifies a scalar curvature singular- 
ity, which is seen for example by inspecting the Kretsch- 
mann invariant R a /s 7 sR a ^ lS ~ e~ 2A ~ exp(2K (i> + uj) 
as v — ► oo. 

Assuming that a similar behaviour occurs in the gene- 
ric case, described by the metric (0) , we can expand Ein- 
stein's vacuum equations R a p = in powers of exp(— A). 
We can then try to solve the zeroth order equations and 
hope that this solution indeed satisfies the assumptions 
made in the expansion. 

Performing such an expansion for the metric (|l|) yields 
equations in u and v only for the metric functions [ [To| , 
i.e. all (^-derivatives are multiplied by at least a factor 
of exp(A). Furthermore it turns out that for Price law 
initial conditions the functions s a ~ exp(A) [ftO| so that 
it effectively drops out of the zeroth order equations. 
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Thus to a very good approximation the vicinity of the 
Cauchy horizon is described by the metric 



so that 



ds 2 



~2e- x dudv + g ab d9 a dO b . 



(3) 



Now the functions A and g a b only depend on the coordi- 
nates u and v. The metric (^) describes the interaction 
region of two plane gravitational waves. This is not an 
accident, but has a deeper physical meaning. Remember, 
that the Cauchy horizon lies at a finite radius r = r_. 
Thus on length scales I < r_ the inner horizons looks 
"flat" and the gravitational waves can be approximated 
by plane waves. Mathematically this corresponds to re- 
placing r 2 dVL 2 by g a b{u 1 v)dO a d6 h . This fact is also re- 
flected in the way the curvature diverges at the Cauchy 
horizon. Rather than 1/r — > oo we have have the metric 
function exp(— A) — > oo as v — > oo. But the metric (^) 
is only a valid approximation if the two-metric g a b stays 
regular (in a sense defined below). Thus for the analysis 
of the full problem (described by the metric (|l|)) to be 
self-consistent the simple plane wave model must show 
the mass inflation behaviour. 

We are therefore left with a two dimensional charac- 
teristic initial value problem in u and v with initial data 
placed on the characteristics S 11 : u = uq and S u : v = Vq 
(see Figure 1). 

Einstein's equations for the metric (H) are pi[ : 

= e A R ab = (4) 
4cr„ d ( a o'„^ - 2d( u (j v ) ab + g ab (K U K V + d( u K v )) 

= R U v = d(uKy) ~ \uv ~ < 7 uab< 7 v ab ~ „K U K V (5) 



and 



= R VV = K v d v \-a vab a v ab --K 2 ~d v K v , (6) 
= Run = K u d u \ - a uab a u ab - X -K 2 U - d u K u (7) 

where lower case indices are raised and lowered with the 
two- metric g ab . The trivial equation (|^) follows from 
equations (||) or (Q) by simple differentiation. 
The quantity 



is the shear of the congruence £^dx a = du A , 



K A = d A yfg 



(8) 



(9) 



is its expansion and we have defined g :— det(g ab ). These 
quantities are again discussed in more detail in [ p"l"| . 

Taking the trace g ab R a b of equation (^) gives d A K A + 
K A K A = or, using the definition (^) for K A , 



P-=V9 = Po + F(v)+G(u). 



(11) 



Comparing this expression to the behaviour of r in 
the spherical case (0) we set 



Po 



-9+1 



-p+1 



(12) 



Note that other choices of the functions F and G merely 
correspond to different gauge choices of the null coordi- 
nates u and v. 

Once this choice has been made one can specify the 
boundary data for the shear a A on the characteristics Yi A 
which comprises the free gravitational data. 

The traceless part of equation (0) is 



^ <J ud{a (J v d 



b) 



0. 



(13) 



Once we have solved this equation for a Aa b we can 
use the constraint equation (||) or (0) to solve for A. 

To be more specific let us choose a definite form for 
the two-metric g a b- A convenient parameterization is fl2]l 



e 2 @ cosh(7) sinh(7) 
sinh(7) e _2/3 cosh(7) 



Sab = P 



The matrix equation ( |13[) then is equivalent to the 
two coupled nonlinear PDE's 



A. 



i 



{P,uP,v + P,vP,u) 



l.uv + — {l.uP.v + l.vP,u) 

2p 



tanh(7) (f3. u j. v + P.vJ.u) 
(14) 

= 2srnh(27)/3 u /3,. y . 



The constraint equations (^|) and (Q) are ordinary first 
order differential equations 



A.4 = 



7J {p> aa 



(p.aY 

2p 



+ 2pcosh 2 7 03 A ) 2 + f(7 A ) 2 ) 



for A — u,v. It is convenient to define a new function 
A := A + 1/2 In p — p tV — p u in terms of which the above 
equations become 



A a = — f2cosh 2 7 (/3. j4 ) 2 + \{ ltA )' 



(15) 



Note that for our choice of p the right hand side of equa- 
tion ( |l5| ) is always smaller or equal to zero for A = v. 
Thus we already can conclude that the function —A grows 
without bound in the limit in oo as expected from the 
spherical case. 

Let us now discuss the boundary data for (3 and 7 
which carry the gravitational degrees of freedom. If we 
ignore the nonlinear terms in the wave equations for [3 
and 7, equations (^) can be rewritten as 



dudvy/g = 



(10) 



□/3 = n 7 = 



3 



where □ is the wave operator corresponding to the metric 
(||). Thus Price's analysis || tells us that the initial data 
on S" is 

a«~7,«~ i-«r p/2 

and on E" is 

where as before q — p + 2. 

Inserting this into equations ( |l5| ) and solving for an 
asymptotic series indeed gives, at least on the initial sur- 
faces 

A „ = + 0(l/v) and A <u = -k\ + 0(l/u) 

where and k? u are constants. This is the expected 
behaviour for A. 

This behaviour was indeed found before by Brady and 
Chambers in jl3| where Einstein's equations for the full 
metric (|l|) were integrated on a pair of intersecting null 
surfaces. 

The important question however is if this behaviour 
still occurs away from the initial surfaces. In our simple 
model this is solely determined by the behaviour of the 
functions (3 and 7, away from the initial surfaces. 

In reference [[Io| perturbative arguments are given 
which show that for Price type initial data the shear must 
decay like a Price tail all the way down to the inner hori- 
zon. The proof however works only for Price power law 
initial data. It would be desirable to see how stable the 
mass inflation scenario is towards changes in the initial 
data. A numerical integration of equations (|lj) and (|l5| ) 
gives the possibility to study the behaviour of the func- 
tions [3 and 7 in detail for different types of boundary 
data, e.g Price power law or a Gaussian pulse. It fur- 
thermore provides a direct verification of the pertubative 
agruments presented in reference ]10| . 

3 Numerical solution of the shear 
equations 

The most natural and obvious approach to the problem 
is to numerically integrate the characteristic initial value 
problem for equations (^) and then use the obtained val- 
ues for the functions (3 and 7 to calculate A using (|l5|). 
While this approach has some problems (see below) it 
is nevertheless the most straightforward way to numer- 
ically integrate these equations. The numerical proce- 
dure, which is schematically displayed in the following 
flow chart, is second order in v and first order in u. 

During our numerical experiments we have used Price 
power law initial data for (3 and a power law or Gaussian 
pulse for 7 on the characteristic initial surfaces E" and 
Ti" . We performed runs for various different values of the 
parameters p and q as well as different amplitudes for 



the initial tails. While a Gaussian pulse as initial data 
changes the simple behaviour of A ~ exp(— k^v — K^u) 
it still leads to the qualitatively same result, namely a 
divergence of A as described above. Power law initial 
data for both functions f3 and 7 on the other hand show 
the behaviour expected from the analytical arguments. 



Calculate (3 and 7 on E Ui from f3 >v 
and 7 it , using a second order method. 






From this calculate A on T, Ui using 
the R vv = equation. 






Insert this into R uu to determine the 
numerical accuracy. 






Use the wave equations fll4|) to cal- 
culate (P v ),u and (7„), u . 






Calculate /3, v and 7,„ on Y< u '+ 5u = 
E" i+1 using a first order method. 



i 



This flow chart shows the principal workings of 
the program used to solve equations (|l4|). E" 4 
denotes the surface of constant u = uo + iSu, 
where Su is the "time" step of the numerical grid. 

In Figure 2 below we show A for such an experiment. 

The problem with the characteristic initial value for- 
mulation is the assumption of an already existing portion 
of the Cauchy horizon. We thus excluded by hand the 
possibility of a spacclikc singularity occurring before a in- 
ner horizon could form. This is acceptable insofar as one 
is looking for a self-consistent description of the struc- 
ture of the Cauchy horizon. If there were a spacelike 
singularity the approximation leading from the general 
metric (|l|) to the plane wave metric (||) would be invalid. 
On the other hand it would be satisfying to know that 
equations ([li]) do not lead to a spacelike singularity for 
generic Price power law initial conditions, placed on a 
spacelike initial surface. 
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To numerically check this one has to change from the 
null coordinates u and v to a set of space- and timelike 
coordinates. This is actually quite easy in the present sit- 
uation thanks to equation ( |io| ) for p. This equation is the 
integrability condition for the coordinate transformation 



P 

X : 



Po 



F(v) 
F(v) 



G(u) 
G(u) 



where F and G have been specified in (12) 
forms the line element (ft) into 



This trans- 



ds z = 



VP 



(d X 2 - dp 2 ) + g ab d9 a d9 b 



(17) 



Note that this coordinate transformation only works be- 
cause p = y/g is harmonic. In other words we can only 
perform this transformation near the Cauchy horizon, 
where the metric (||) is valid. Thus, even though the 
transformation is defined for v > and — u > it only 
is applicable near the Cauchy horizon, i.e. in the limit of 
large v. 




Figure 2. Shown is the function A which corre- 
sponds to the negative logarithm of the "mass 
function" of the hole. Thus m increases expo- 
nentially towards the inner horizon (v —* oo). 

Initial data is now placed on the surface p = po- (Note 
that p, which in the spherical case is proportional to the 
radial coordinate r 2 is timelike, see Figure 3.) Only this 
surface goes into the "corner" H of Figure 1. For p < po 
it intersects the inner horizon while, for p > po it seems 
to intersect the event horizon. However the plane wave 
metric (^|) is not a valid approximation of the spacetime 
in the limit u — > —oo, v = constant. 

The equations governing the functions 7 and (3 in 
terms of the new coordinates are 



P,pp P,xx + 



P 



2tanh( 7 ) (P, x % x - 0, P %, 



7,pp - 7,xx + — = 2 si!1 



h(2 7 ) ((/3 P ) 2 - {pj 



The two constraint equations for A now become 

2pcosh 2 ( 7 )((/? p ) : 
+ 1 ((7,p) 



(7, x ) 2 ) 



(18) 



(16) and 



A 



4pcosh 2 (7)/? x /? p + (n, x l,, 




Figure 3. This Figure shows the new coordinates 
p and x- Initial data is now placed on the surface 
p — po, and the evolve to smaller values of p, i.e. 
towards the inner horizon. 



Again, because the right hand side of equation (|18|) is 
always greater or equal to zero A has to decrease as we 
go to smaller values of p, i.e. towards the Cauchy hori- 
zon. While we still expect the functions f3 and 7 to stay 
bounded this is not true any more for their derivatives, 
since terms like dv/dp which come into the coordinate 
transformation blow up. Thus we still expect A to di- 
verge as we approach the inner horizon. 

To integrate these equations we have used a centered 
two step Lax-Wendrof scheme H] which is second order 
in p and x- The program now evolves the initial data 
inside the domain of dependance of the initial Cauchy 
surface, so that no assumptions about boundary data 
have to be made. For appropriate values of the numerical 
grid constants this domain covers almost the entire region 
between the surface p — po and the inner horizon. 

Again the results agree with the analytical approx- 
imations as well as the characteristic numerical results 
indicating the validity of the characteristic code. In par- 
ticular we don't see any spacelike singularity developing. 
It must be added however that even this approach is not 
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completely useful in determining this question. While the 
initial surface p = po approaches the "corner" H it also 
approaches the inner horizon, which initially is at radius 
p = Pq. Ideally one would like to deal with a space like 
Cauchy surface that never touches the inner horizon, but 
still always stays below the potential barricrc. However 
this cannot be accomplished by the simple coordinate 
transformation (|l6|). On the other hand the two surfaces 
only meet at infinitly late time u — — oo and v = oo. 

4 Conclusion 

The main purpose of this paper is to support and sup- 
plement the analytical work done |, §, on the str uc- 
ture of the Cauchy horizon or mass inflation singularity. 
While the latter is necessarily confined to perturbative 
methods a numerical analysis gives one the opportunity 
to investigate specific examples for nontrivial realistic ini- 
tial data. 

The results of our analysis strongly support the result 
of the analytical work. We indeed find the predicted 
behaviour of the shear and the metric function A, which 
is essentially the negative logarithm of the (quasi-local 
15 1) mass function of the black hole. It thus seems that 
mass inflation singularity is a generic singularity [jlo| . [k]] . 
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